run scripts/start_RioTinto.m

%% Import data
opts = spreadsheetImportOptions("NumVariables", 5);
addpath('figures')

% Specify sheet and range
opts.Sheet = "ForMatlabExport";
opts.DataRange = "A2:E1512";

% Specify column names and types
opts.VariableNames = ["Dates", "PriceLtd", "VolumeLtd", "PricePLC", "VolumePLC"];
opts.VariableTypes = ["string", "double", "double", "double", "double"];

% Specify variable properties
opts = setvaropts(opts, "Dates", "WhitespaceRule", "preserve");
opts = setvaropts(opts, "Dates", "EmptyFieldRule", "auto");

% Import the data
tbl = readtable("raw_data\RioTintoData.xlsx", opts, "UseExcel", false);

%% Convert to output type
Dates     = tbl.Dates;
PriceLtd  = tbl.PriceLtd;
VolumeLtd = tbl.VolumeLtd;
PricePLC  = tbl.PricePLC;
VolumePLC = tbl.VolumePLC;

clear opts tbl

DateVec = Dates(1:end,:);
DateVec=datetime(DateVec);
PriceLtd = PriceLtd(1:end,:);
VolumeLtd = VolumeLtd(1:end,:);
PricePLC = PricePLC(1:end,:);
VolumePLC = VolumePLC(1:end,:);

meanDVs = zeros(size(Dates,1),2);
counter = 0;
while (counter<size(Dates,1))
   counter = counter + 1; 
   if(counter<261)
       fracDVs(counter,:) = [0 0];
   else
       %Multiply by 10 because of scaling: prices are rescaled up by 100
       %and volumes are scaled down by 1000 so 1000/100 = 10 for true value
       fracDVs(counter,:) = [mean(VolumeLtd(counter-260:counter,1).*(PriceLtd(counter-260:counter,1))) mean(VolumePLC(counter-260:counter,1).*(PricePLC(counter-260:counter,1)))];
   end
end

arbSizes = zeros(size(Dates,1)-260,1);
totalWelfare = zeros(size(Dates,1)-260,1);
fnEvals = zeros(size(Dates,1)-260,1);

counter = 260;
while(counter<size(Dates,1))
   counter = counter + 1;
   Values = @(m) (PriceLtd(counter,1).*(1+(100*(1000000*m/fracDVs(counter,1))*0.0000) + (((100*abs(1000000*m)./fracDVs(counter,1)).^0.5)*0.000889.*sign(m))))-(PricePLC(counter,1)*(1-(100*(1000000*m/fracDVs(counter,2))*0.0000) - (((100*abs(1000000*m)/fracDVs(counter,2)).^0.5)*0.000889.*sign(m))));
   Values2 = @(m) (PriceLtd(counter,1).*(1+(100*(1000000*m/fracDVs(counter,1))*0.0000) + (((100*abs(1000000*m)./fracDVs(counter,1)).^0.5)*0.000889.*sign(m))))./(PricePLC(counter,1)*(1-(100*(1000000*m/fracDVs(counter,2))*0.0000) - (((100*abs(1000000*m)/fracDVs(counter,2)).^0.5)*0.000889.*sign(m))))-1;
   
   options = optimoptions('fsolve','FunctionTolerance',1e-12,'OptimalityTolerance',1e-12,'StepTolerance',1e-6,'MaxFunctionEvaluations',10000);
   [m1,fnVal] = fsolve(Values,1,options);
   grid = sign(m1)*(0:1:abs(m1));
   evals2 = zeros(size(grid,2),1);
   counter2 = 0;
   while(counter2<size(grid,2))
       counter2 = counter2+1;
       evals2(counter2,1) = Values2(grid(1,counter2));
   end
   arbSizes(counter-260,1) = m1*1000000;
   %evals = Values(grid);
   totalWelfare(counter-260,1) = sum(evals2)*1000000;
   fnEvals(counter-260,1) = fnVal;
end

ratios = PriceLtd./(PricePLC);
ratios = ratios(261:end,:);

save('intermediate\workspace_RioTinto');

%% Create figures
load('intermediate\workspace_RioTinto');

fig_RioTinto_welfare(DateVec, ratios, totalWelfare, optfig)
fig_RioTinto_arbsize(DateVec, ratios, arbSizes, optfig)
